﻿module PCA

open System
open System.IO
open System.Collections.Generic
open MathNet.Numerics
open MathNet.Numerics.FSharp
open System.Diagnostics
open MathNet.Numerics.LinearAlgebra.Double
open Microsoft.Office.Interop.Excel

//set point as default decimal separator
System.Threading.Thread.CurrentThread.CurrentCulture <- System.Globalization.CultureInfo.CreateSpecificCulture("En")
System.Threading.Thread.CurrentThread.CurrentUICulture <- System.Globalization.CultureInfo.CreateSpecificCulture("En")

//let eigenvectors (matrix:DenseMatrix) =
//  let eigen = matrix.Evd()
//  let evals = eigen.EigenValues()
//  let res = DenseMatrix(matrix.RowCount,matrix.ColumnCount)
//  for r = 1 to matrix.RowCount do
//    //Console.WriteLine(evals.[r - 1].Real * a.Row(r - 1))
//    res.SetRow(r - 1,evals.[r - 1].Real * matrix.Row(r - 1))
//  res

type ImageData = 
  {
    X : DenseVector
    Y : DenseVector
  }

let text = File.ReadAllText("semeion.data")
let dimensions = 256
let mutable optimum = false
Console.WriteLine("Insert the number of times the knn is run to compute the mean error")
let test_number = Console.ReadLine() |> Convert.ToInt32
Console.WriteLine("Would you like to build an output excel data sheet? (Enter 0, otherwise a text file will be created instead)")
let mutable excel_output = true
let mutable output_directory = ""
if Console.ReadLine() <> "0" then
  excel_output <- false
  Console.WriteLine("Insert the directory of the output file")
  output_directory <- Console.ReadLine()
  if output_directory.EndsWith("\\") |> not then
    output_directory <- output_directory + "\\"

Console.WriteLine("Initializing data... Please wait...")
let dataList = text.Split(' ','\n') |> Seq.toList |> List.filter(fun s -> s <> "" && s <> "\r")
let (data_vector:List<ImageData>) = ResizeArray()
let seed = Random()

//gather data and put them in a list of ImageData
let rec gather_data (data:string list) counter x_vector =
  let rec gather_labels (l: 'a list) steps res =
    if steps = 0 then
      (l,res)
    else
      gather_labels l.Tail (steps - 1) (res @ [l.Head])

  if data = [] then
    ()
  else
    if counter = dimensions then
      let l,y = gather_labels data 10 []
      let labels = y |> List.map(fun y1 -> y1 |> Double.Parse)
      let image_component =
        {
          X = DenseVector(x_vector |> List.toArray)
          Y = DenseVector(labels |> List.toArray)
        }
      data_vector.Add(image_component)
      gather_data l 0 []
    else
      gather_data data.Tail (counter + 1) (x_vector @ [data.Head |> Convert.ToDouble])

gather_data dataList 0 []
let n = data_vector |> Seq.length


//compute the mean vector
let xm (data_vector:List<ImageData>)=
  let mutable sum = DenseVector(dimensions)
  for i = 1 to data_vector.Count do
    sum <- sum + data_vector.[i - 1].X
  sum / (float data_vector.Count)

//compute the covariance matrix
let cov_matrix (xm:DenseVector) data_vector =
  let n = data_vector |> Seq.length
  let mutable res = DenseMatrix(dimensions)
  for image in data_vector do
    res <- res + (DenseVector.OuterProduct((image.X - xm),(image.X - xm)))
  res <- res.Divide(float n):?>DenseMatrix
  res

//compute the eigenvectors
let eigen_vectors (matrix:DenseMatrix) =
  let eigen = matrix.Svd(true)
  let ev = eigen.VT():?>DenseMatrix
//  let sorted_ev = DenseMatrix(ev.RowCount,ev.ColumnCount)
//  //reverse the matrix since the eigenvalues are in increasing order and paired with the corresponding eigenvectors
//  for i = 0 to ev.RowCount - 1 do
//    sorted_ev.SetRow(i,ev.Row(ev.RowCount - 1 - i))
//  sorted_ev
  ev

//compute the principal components vector
let alpha (vectors:ResizeArray<ImageData>) (eigen_vectors:DenseMatrix) xm =
  let alpha_mat = DenseMatrix(vectors.Count,eigen_vectors.RowCount)
  for i = 1 to vectors.Count do
    for j = 1 to eigen_vectors.RowCount do
      alpha_mat.[i - 1,j - 1] <- (vectors.[i - 1].X - xm) * eigen_vectors.Row(j - 1)
  alpha_mat

//assign a label to each class
let label_id (x:DenseVector) =
  let mutable id = 0
  for i = 1 to 10 do
    if x.[i - 1] = 1.0 then 
      id <- i
  id

//create a list containing a vector of points divided into classes
let feature_split_data = ResizeArray()
for i = 1 to 10 do
  let split_list =
    List(seq{for image in data_vector do
              if label_id image.Y = i then yield image})
  feature_split_data.Add(split_list)
  

//function to pick a random element from a list
let pick_random_subset (list:List<'a>) (res:List<'a>) size =
  for i = 1 to size do
    let picked_element = list.[seed.Next(0, list.Count)]
    list.Remove(picked_element) |> ignore
    res.Add(picked_element)

//function picking 10 random elements per class
let build_data (training_set:List<ImageData>)=
  for i = 1 to 10 do
    let rand_set = ResizeArray()
    for x in feature_split_data.[i - 1] do
      rand_set.Add(x)
    pick_random_subset rand_set training_set 10

//Principal Component Analysis function
let pca (data_set:seq<ImageData>) training_set xm_data xm_train =
  let ev = eigen_vectors (cov_matrix xm_train training_set)
  let alpha_training = alpha training_set ev xm_train
  let alpha_data = alpha (ResizeArray(data_set)) ev xm_train
  (alpha_training,alpha_data)

let ones n =
  let out = DenseMatrix(n,1)
  for i = 1 to n do
    out.[i - 1,0] <- 1.0
  out
  
let sum_column (a:DenseMatrix) =
  let out = DenseMatrix(1,a.ColumnCount)
  for j = 1 to a.ColumnCount do
    let mutable sum = 0.0
    for i = 1 to a.RowCount do
      sum <- sum + a.[i - 1,j - 1]
    out.[0,j - 1] <- sum
  out

let matrix_square (a:DenseMatrix) =
  let out = DenseMatrix(a.RowCount,a.ColumnCount)
  for i = 1 to a.RowCount do
    for j = 1 to a.ColumnCount do
      out.[i - 1,j - 1] <- a.[i - 1,j - 1] * a.[i - 1, j - 1]
  out

(*function computing the euclidean distance matrix between two matrices: if the first matrix
X is a AxB matrix and Y is a CxD matrix then the output matrix is a AxC matrix*)
let distance_matrix (a:DenseMatrix) (b:DenseMatrix) =
  let ob = ones(b.RowCount)
  let sa = sum_column ((matrix_square a).Transpose() :?> DenseMatrix)
  let x = (ob * sa).Transpose()
  let y =  ones(a.RowCount) * sum_column((matrix_square b).Transpose() :?> DenseMatrix)
  let z = 2.0 * a * (b.Transpose()) :?> DenseMatrix
  x + y - z :?> DenseMatrix

//function to compute K-nn classification algorithm
let knn components (a_data:DenseMatrix) (a_training:DenseMatrix) (training_set:List<ImageData>) =
  let reduced_training = a_training.SubMatrix(0,a_training.RowCount,0,components) :?> DenseMatrix
  let reduced_data = a_data.SubMatrix(0,a_data.RowCount,0,components) :?> DenseMatrix
  let dist = distance_matrix reduced_data reduced_training
  seq{for i = 1 to data_vector.Count do
        let min = dist.Row(i - 1).MinimumIndex()
        yield training_set.[min].Y}

//function which tests knn for a number of components which varies from 1 to 256
let test_knn a_data a_training data_set (training_set:List<ImageData>) =
  let error_vector = DenseVector(256)
  for i = 1 to 256 do
    let mutable total_error = 0
    let result = knn i a_data a_training training_set
    let result_list = Seq.zip data_set result
    for (data,classification) in result_list do
      if data.Y <> classification then
        total_error <- total_error + 1
    error_vector.[i - 1] <- (float total_error) / (data_set |> Seq.length |> float)
  error_vector

Console.WriteLine("Done!")
Console.WriteLine("Computing KNN... Please wait...")
let timer = new Stopwatch()
let error_matrix = DenseMatrix(test_number,256)

//*** EXCEL DATA TABLE ***
let build_excel_table (a_data:DenseMatrix) (mean_error_vector:DenseVector) =
  // Run Excel as a visible application
  let app = new ApplicationClass(Visible = true) 

  // Create new file and get the first worksheet
  let workbook = app.Workbooks.Add(XlWBATemplate.xlWBATWorksheet) 
  // Note that worksheets are indexed from one instead of zero
  let worksheet = (workbook.Worksheets.[1] :?> Worksheet)

  // Store data in arrays of strings or floats
  let rnd = new Random()
  let error_data = Array2D.init mean_error_vector.Count 2 (fun i j ->
                                                                if j = 0 then (i + 1) |> float else (1.0 - mean_error_vector.[i]))
  let component_data = Array2D.init a_data.RowCount 2 (fun i j -> a_data.[i,j])
  let component_class = Array2D.init data_vector.Count 1 (fun i _ -> (label_id data_vector.[i].Y))

  // Populate data into Excel worksheet
  worksheet.Range("A1", "B256").Value2 <- error_data
  worksheet.Range("C1","D" + (a_data.RowCount).ToString()).Value2 <- component_data
  worksheet.Range("E1","E" + (data_vector.Count).ToString()).Value2 <- component_class

let build_output_file (a_data:DenseMatrix) (mean_error_vector:DenseVector) =
  Directory.CreateDirectory output_directory |> ignore
  File.WriteAllLines(output_directory + "pca.dat",["Alpha first two components:\n\n"] @
    [a_data.ToString()] @
    ["\n\nError values for each class:\n\n"] @
    [mean_error_vector.ToString()])

//build a matrix of error vectors
timer.Start()
let mutable alpha_data_matrix = DenseMatrix(data_vector.Count,256)
for i = 1 to error_matrix.RowCount do
  let training_set = ResizeArray()
  build_data training_set
  let data_set = data_vector |> Seq.filter(fun x -> training_set |> Seq.exists(fun x1 -> x1 <> x))
  let xm_train = xm training_set
  let xm_data = xm (List(data_set))
  let (a_training,a_data) = pca data_set training_set xm_data xm_train
  alpha_data_matrix <- a_data
  error_matrix.SetRow(i - 1,test_knn a_data a_training data_set training_set) //each error matrix contains the error vector for each instance of knn
  Console.WriteLine("Data set: " + i.ToString() + " succesfully computed!")
  Console.WriteLine("Done!")

//compute the mean component of the vectors in a matrix
let mean_vector (matrix:DenseMatrix) =
  let mean_vector = DenseVector(matrix.ColumnCount)
  for j = 1 to matrix.ColumnCount do
    let currentColumn = error_matrix.Column(j - 1)
    let mutable sum = 0.0
    for i = 1 to currentColumn.Count do
      sum <- sum + currentColumn.[i - 1]
    mean_vector.[j - 1] <- sum / (float error_matrix.RowCount)
  mean_vector

Console.WriteLine("Generating output file...")
if excel_output then
  build_excel_table alpha_data_matrix (mean_vector error_matrix)
else
  build_output_file alpha_data_matrix (mean_vector error_matrix)
Console.WriteLine("Done!")


timer.Stop()
Console.WriteLine("Elapsed time: " + timer.ElapsedMilliseconds.ToString())


